from eo_utils import *
Before loading and processing the data, we can have a look at the available collections:
backend = "https://openeo.cloud"
con = openeo.connect(backend)
con.list_collections()
Get the metadata for Sentinel-2 ARD collection:
con.describe_collection("boa_sentinel_2")
Authenticate using EGI
conn = openeo.connect(backend).authenticate_oidc(provider_id='egi')
center = [47.29, 10.75]
zoom = 7
eoMap = openeoMap(center,zoom)
addS2Tiles(eoMap)
eoMap.map
Get the bounding box from the previous map.
bbox = eoMap.getBbox()
print('west',bbox[0],'\neast',bbox[2],'\nsouth',bbox[1],'\nnorth',bbox[3])
west 11.066666 east 11.171036 south 46.038835 north 46.110286

collection = 'TERRASCOPE_S2_TOC_V2'
spatial_extent = {'west':bbox[0],'east':bbox[2],'south':bbox[1],'north':bbox[3],'crs':'EPSG:4326'}
temporal_extent = ["2020-08-08", "2020-08-12"]
bands = ["B04"]
s2_l2a = conn.load_collection(collection,spatial_extent=spatial_extent,bands=bands,temporal_extent=temporal_extent)
s2_l2a_geotiff = s2_l2a.reduce_dimension(reducer="mean", dimension="t").save_result(format="GTIFF")
job = conn.create_job(s2_l2a_geotiff)
job_id = job.job_id
if job_id:
print("Batch job created with id: ",job_id)
job.start_job()
else:
print("Error! Job ID is None")
Batch job created with id: vito-66cc8836-1370-4682-8332-10a8ce9bbc21
We can get a description of the job and check its status.
job = conn.job(job_id)
job_description = job.describe_job()
print("Batch job with id: ",job_id, ' is ',job_description['status'])
Batch job with id: vito-66cc8836-1370-4682-8332-10a8ce9bbc21 is running
job
results = job.get_results()
results
results.download_file("./data/US1/s2_result.tif")
Running the next cell, the Sentinel-2 data will appear on the map as in the following screenshot.

addLayer(eoMap,'./data/US1/s2_result.tif','S2 B04',clip=[0,1500])
eoMap.map
B04_S2 = xr.open_rasterio('./data/US1/s2_result.tif')
B04_S2
<xarray.DataArray (band: 1, y: 829, x: 965)>
[799985 values with dtype=uint16]
Coordinates:
* band (band) int64 1
* y (y) float64 5.156e+06 5.156e+06 5.156e+06 ... 5.148e+06 5.148e+06
* x (x) float64 6.743e+05 6.743e+05 6.743e+05 ... 6.839e+05 6.839e+05
Attributes:
transform: (10.0, 0.0, 674250.0, 0.0, -10.0, 5156280.0)
crs: +init=epsg:32632
res: (10.0, 10.0)
is_tiled: 1
nodatavals: (0.0,)
scales: (1.0,)
offsets: (0.0,)
AREA_OR_POINT: Area[799985 values with dtype=uint16]
array([1])
array([5156275., 5156265., 5156255., ..., 5148015., 5148005., 5147995.])
array([674255., 674265., 674275., ..., 683875., 683885., 683895.])
fig, ax1 = plt.subplots(1,1,figsize=(20,20))
ax1.imshow(B04_S2[0].values,cmap='Greys_r',vmin=0,vmax=2000)
<matplotlib.image.AxesImage at 0x7fb4911ee160>
collection = 'TERRASCOPE_S1_GAMMA0_V1'
spatial_extent = {'west':bbox[0],'east':bbox[2],'south':bbox[1],'north':bbox[3],'crs':'EPSG:4326'}
temporal_extent = ["2019-09-18", "2019-09-25"]
bands = ["VH"]
s1 = conn.load_collection(collection,spatial_extent=spatial_extent,bands=bands,temporal_extent=temporal_extent)
s1_geotiff = s1.reduce_dimension(reducer="max", dimension="t").save_result(format="GTIFF")
s1_job = conn.create_job(s1_geotiff)
s1_job_id = s1_job.job_id
if s1_job_id:
print("Batch job created with id: ",s1_job_id)
s1_job.start_job()
else:
print("Error! Job ID is None")
Batch job created with id: vito-9879f280-2f9d-40b1-950c-81a5263296a5
s1_job = conn.job(s1_job_id)
job_description = s1_job.describe_job()
print("Batch job with id: ",s1_job_id, ' is ',job_description['status'])
Batch job with id: vito-9879f280-2f9d-40b1-950c-81a5263296a5 is running
results = s1_job.get_results()
metadata = results.get_metadata()
results.download_file("./data/US1/s1_result.tif")
WindowsPath('data/US1/s1_result.tif')
Running the next cell, the Sentinel-1 data will appear on the map as in the following screenshot.

addLayer(eoMap,'./data/US1/s1_result.tif','VH',clip=[0,0.1])
eoMap.map
We can have a look at the original spatial coverage of the ARD Landsat-8 collection, with its four overlapping areas.
We are going to demonstrate that, with the Equi7 grid used in this ARD collection, the data is being reprojected and therefore there are no problems of duplicated data in overlapping areas.
conn_2 = openeo.connect("https://openeo.eodc.eu/v1.0")
itemIter = conn_2.collection_items(
"boa_landsat_8",
temporal_extent = ["2019-09-20", "2019-09-22"],
spatial_extent = [5.0,45.0,20.0,50.0]
)
next(itemIter)
center = [45.85, 11.0]
zoom = 11
eoMap2 = openeoMap(center,zoom)
eoMap2.map
bbox = eoMap2.getBbox()
print('west',bbox[0],'\neast',bbox[2],'\nsouth',bbox[1],'\nnorth',bbox[3])
west 10.809517 east 10.91526 south 45.850916 north 45.9121
collection = 'boa_landsat_8'
spatial_extent = {'west':bbox[0],'east':bbox[2],'south':bbox[1],'north':bbox[3],'crs':'EPSG:4326'}
temporal_extent = ["2019-09-20", "2019-09-22"]
bands = ["B02","B03","B04"]
l8 = conn.load_collection(collection,spatial_extent=spatial_extent,bands=bands,temporal_extent=temporal_extent)
l8_nc = l8.save_result(format="netCDF")
l8_job = conn.create_job(l8_nc)
l8_job_id = l8_job.job_id
if l8_job_id:
print("Batch job created with id: ",l8_job_id)
l8_job.start_job()
else:
print("Error! Job ID is None")
l8_job = conn.job(l8_job_id)
job_description = l8_job.describe_job()
print("Batch job with id: ",l8_job_id, ' is ',job_description['status'])
results = l8_job.get_results()
metadata = results.get_metadata()
results.download_file("./data/US1/l8_result2.nc")
PosixPath('data/US1/l8_result2.nc')
dc = xr.open_dataset('./data/US1/l8_result2.nc')
dc
<xarray.Dataset>
Dimensions: (y: 274, x: 310)
Coordinates:
time datetime64[ns] 2019-09-21T12:00:00
spatial_ref int32 0
* y (y) float64 1.424e+06 1.424e+06 ... 1.416e+06 1.416e+06
* x (x) float64 4.815e+06 4.815e+06 ... 4.825e+06 4.825e+06
Data variables:
B02 (y, x) float64 190.0 192.0 181.0 193.0 ... 263.0 272.0 264.0
B03 (y, x) float64 369.0 359.0 330.0 373.0 ... 548.0 550.0 541.0
B04 (y, x) float64 226.0 230.0 213.0 235.0 ... 296.0 311.0 296.0
Attributes:
crs: PROJCS["Azimuthal_Equidistant",GEOGCS["WGS 84",DATUM[...
datetime_from_dim: 2019-09-21T12:00:00.000000000array('2019-09-21T12:00:00.000000000', dtype='datetime64[ns]')array(0, dtype=int32)
array([1424445., 1424415., 1424385., ..., 1416315., 1416285., 1416255.])
array([4815375., 4815405., 4815435., ..., 4824585., 4824615., 4824645.])
array([[190., 192., 181., ..., 391., 409., 374.],
[183., 186., 177., ..., 356., 295., 298.],
[181., 176., 172., ..., 424., 388., 371.],
...,
[303., 285., 263., ..., 263., 274., 308.],
[332., 279., 246., ..., 272., 266., 267.],
[441., 281., 264., ..., 263., 272., 264.]])array([[369., 359., 330., ..., 546., 591., 560.],
[349., 350., 332., ..., 550., 475., 497.],
[340., 341., 339., ..., 566., 493., 474.],
...,
[406., 420., 411., ..., 557., 572., 592.],
[503., 436., 376., ..., 561., 553., 562.],
[594., 461., 408., ..., 548., 550., 541.]])array([[226., 230., 213., ..., 505., 563., 547.],
[218., 221., 208., ..., 472., 385., 400.],
[222., 219., 220., ..., 551., 487., 457.],
...,
[299., 322., 299., ..., 320., 345., 410.],
[367., 316., 273., ..., 304., 302., 329.],
[578., 374., 321., ..., 296., 311., 296.]])We can visualize the RGB composite as a simple image using Matplotlib. Notice the seamless integration of the overlapping Landsat tiles.
fig, ax1 = plt.subplots(1,1,figsize=(10,10))
ax1.imshow(tone_mapping(dc.B04,dc.B03,dc.B02),cmap='brg')
<matplotlib.image.AxesImage at 0x7fb492481a90>